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Abstract 

In this work we consider topology optimization of systems, which are governed by the exter- 
nal Bernoulli free boundary problem. We utilize the so-called pseudo-solid approach to solve the 
governing free boundary problems during the optimization. To define design domains we utilize 
a level set representation parameterized by radial basis functions. This design parametrization 
allows topological changes in the design domain. 

Introduction 

In shape optimization, state problems are usually given by partial differential equations which are for- 
mulated and solved in known domains. Thus the unknown shape we are looking for is represented by 
a solution to our optimization problem. If the state is defined by a free boundary problem, the situation 
becomes more involved. The main feature of such problems is the fact that not only a function solv- 
ing the respective PDE but also the domain itself where this PDE is considered are unknown and they 
have to be found simultaneously. Unlike to the standard case, the unknown geometry now appears in 
both, the upper as well as the inner level of the optimization problem. Classical shape optimization 
with external Bernoulli free boundary problems (BFBP) as the state relation has been presented in 
im. A couple {u, Q) is said to be a solution to BFBP if u is harmonic in a doubly connected (un- 
known) domain Q, it satisfies the Dirichlet condition on the boundary dco of the inner component uj 
given a-priori and, in addition the over-determined system consisting of the Neumann and Dirichlet 
condition on the (unknown) free boundary. The shape of the free boundary was controlled by lo but 
still keeping the same topology of co. To make numerical realization simpler we supposed that the 
system of admissible wis consists of star-like domains enabling us to express duj and in terms of 
the polar coordinates. However it was observed that this system is very narrow and many reasonable 
target free boundaries can not be matched. Indeed, if the boundary du is twice differentiable and co is 
star-like with respect to a neighborhood of some point in co then the respective free boundary is of the 
class C°° in. Hence if a given target free boundary has the discontinuous curvature (a square with 
rounded comers, e.g) then it can never be realized for any such uj. In computations this fact manifests 
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itself by oscillations of duo: using finer discretizations of du, the oscillations become faster and faster 
(see Q). This behavior indicates the tendency to change the topology of w. There are two ways how 
to overcome such oscillations: i) to restrict the design space ii) to extend it in such a way that changes 
of topology are possible. If i) is used then oscillations are suppressed but usually there is a big gap 
between the target and the found free boundary. To get better results one has to change the topology 
of ijj. One of possible ways how to do that will be described in this paper. We shall use a level set ap- 
proach. The topology of w will determined by the zero level set of a function whose argument is given 
by the linear combination of a system of radial basis functions with overlapping supports. Instead of 
solving Hamilton- Jacobi equation describing the evolution of the level set function, the problem can 
be treated as a parametric minimization problem with parameters represented by the coefficients of 
the above mentioned linear combination. 

The paper is organized as follows: Section 1 presents a general setting of a class of topology 
optimization problems governed by BFBP. Section 2 is devoted to the presentation of the state solver. 
Since the state problem will be solved several times, one has to have at his disposal an efficient and 
reliable method. It turns out that the so-called pseudo-solid approach enjoys both these requirements. 
In Section 3 we shortly recall a level set approach widely used in topology optimization. Section 
4 deals with a discretization of the whole optimization problem and its numerical realization when 

radial basis functions are used to parameterize the level set function. Section [s] presents results 
of several model examples. Finally, the paper is completed with two appendices on smoothing the 
Heaviside function and on the analytical solution to a specific Bernoulli problem. 

Throughout the paper we use the following notation: the symbol H^{Q) {k > integer) stands 
for the Sobolev space of functions which are together with their derivatives up to order k square 
integrable in Q, i.e. elements of L'^{Q) (we set = L^($7)). 



1 Setting of the problem 



We start with the definition of the state problem represented by an exterior Bernoulli free boundary 
problem. Let 7 < and an open set a; C with a sufficiently regular boundary duj be given. The 
problem consists in finding a set 17 D cJ and a function n : J] \ a; — )• M satisfying 

An = in \ cj 
u = 1 on d(jj 




on dQ. 



This paper deals with the control of the shape of in {Vioj) I. The geometry of the free boundary 
do. will be driven by the shape of vo towards a given target free boundary Tt- Our optimization 
problem reads as follows: 

f Find bj* such that 

< (P) 

\ J(r(a;*)) < J(r(^)) 

holds for any w G O, where O is a set of admissible designs. The cost functional to be cons idered is 
the distance between Tt and the free boundary r(cj) := being the solution of {V{uj) I, i.e. 

JiTiu;)) = p{T{u;),Tt), (1) 

where p is a function characterizing the distance of r(a;) from Tt- 
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In this paper we make the assumption that $7 is a star-like domain. However, in contrast to lUl no 
such assumption is made on u meaning that also the topology of lo may change. 



2 Pseudo-solid formulation of Bernoulli free boundary problems 



Free boundary problems have in common the difficulty that the geometry (here the domain Q) has to 
be determined simultaneously with the solution u of the state problem, which implies that a numerical 
solution has to be done iteratively [3]. Possible solution strategies include trial methods, linearization 
methods (continuous or discrete) H, and shape optimization methods Q. We decided for the so- 
called pseudo-solid approach (PSA) in which the unknown domain is obtained by an appropriate 
deformation of the reference configuration U. This deformation is one of the unknowns in PSA. The 
main advantage of this approach is that there is no need to construct an explicit parametrization of the 
geometry using e.g. a conformal mapping of the reference domain. PSA is useful also in numerical 
realization: finite element partitions of Q can be constructed via the respective deformation of the 
partition of 17. As the Bernoulli problem has to be solved several times for different uj, the choice 
of an efficient solver is also important. We use Newton's method because of its fast convergence. 
Moreover, this solution strategy readily allows us to obtain geometrical sensitivities of the system. 

Let C be a fixed, simply connected reference domain. In the pse udo-so lid technique we 
construct a mapping F : — )• M^, Q, i— )• F{Q,) =: Q such that Q, solves {V{uj) I for given uj. To 



construct such F we treat Q, as an elastic solid that undergoes a deform ation ca used by an external 
loading p such that the deformed solid defines such Q. Thus, problem {V{lo) i is strongly coupled 
with the Lame system of linear elasticity in which the loading p applied to dQ plays the role of an 
unknown in PSA. This approach has been previously used to solve free surface flow problems (see 
e.g. [6, 7]) and Bernoulli free boundary problems in [1]. 

For any w G VF"*^ :={"sufficiently" small and regular deformations} we define the domain 



{x e 



X 



X + w{x), X G rj}. 



We introduce the following function spaces: 

W^ = {we [H\h)f \ w\^ = 0} 

Vcin) = WeH\n)\ip\^ = c}, cgm. 

The pseudo-solid formulation of our free boundary problem then reads as follows: Given uj, 
find {u,p,v) e Vi{Q.y) X L'^{dn) x W"^ such that 



Vu-Vipdx = j / ipds \/ip£Vo{^v) 



a(v) : e(w) dx 



n\uj 



pn - w ds \/w G . 



an 



Equations (P(a;)i) and {V{uj)2) constitute the weak form of (T'{uj) I while (P(a;)3) is the weak 



form of the linear elasticity problem in ^}\uj. Here pn is an (unknown) external load. The components 
of the strain and stress tensors e = {eij{v)}md a = {(^ij{v)} associated with a displacement field v 
are given by 



1 / dv 



dxj 



hj = 1,2, 
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respectively, where and A are the Lame coefficients. Since in this case the linear elasticity system 
does not have any physical meaning, the Lame coefficients can be chosen quite freely. In this paper 
the choice = 0.5 a nd A = was made. The solvability of I is analyzed in details in HI. 



The relation between (V{ui) I and ^ 
couple (u| , solves (V{u!) 



is readily seen: if {u,p, v) is a solution of {V{uj) I then the 



External load p 




u = and grad(u).n = y 



a) reference configuration b) deformed configuration 

Figure 1 : Principle of the pseudo-solid approach. 



3 Level set approach to the optimization problem 

In the previous paper 111 the design domain was parameterized using the polar co-ordinates. A 
tendency towards fractal-like designs indicating possible topological changes of uj was observed in 
certain cases. However, the boundary variation technique used in [ 1 ] is not able to handle topological 
changes automatically. This will be done by a level set parametrization ||8l|9l of w. 

The basic idea is simple: Let D be a larger domain containing all admissible io (for example a 
rectangle slightly larger than the bounding box of the target boundary Ft). Let ip : D ^ Rbe given 
and define the set u> by 

UJ := uj{ij) :={x£D\ tPix) > 0} , tp £ U""^. (2) 

Here U"''^ is a family of admissible level set functions such that ^ G 11°-'^ =^ '^(V') G O. Clearly 
the parametrization (|2]) allows topological changes of co. 
We can now reformulate problem (P) as follows: 

Find ^b* e W"^ such that 

j(r(^*)) < J{r{ij)) w G u"^, 

where T{Tp) is the free boundary defined by {V{uj{ip))) 

Next we introduce the following relaxed state problem which does not contain explicitly any 
Dirichlet type boundary conditions. For e > 0, G U"''^ given: 
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Find {ue,Pe,v,) E iJ^^^^J x L'^id^) x [H^h)]'^ sucli that 

/ Vue-Vifdx-'y Lpds+ I Ge{ij){ue - l)ip dx = V(/J G 

la{vs):e{w)dx- l peU- wds + lGe{ip)v£-wdx = Vw G 
Jn Jan. Jn 

Above, Ge is a penalty function releasing the constraints n = and i> = in clj. The classical choice 
is Ge = \H, where H is the Heaviside function. This choice of G^ will be used in what follows. 
We define now the "relaxed" optimization problem 

f Find lb* e U"'^ such that 

< (P ) 

I j(r,(v:) < J(r,W) vv^GC/'^^ 



where Teiip) = Fs{dO,), = id + Vs and Vs is the third component of the solution of ('Pe(V') 



4 Discretization 

One of the advantages of level set methods is that they avoid tracking of the boundary of the design 
domain. Instead a fixed mesh is used. In our case the mesh T of the reference domain is fixed, 
but the mesh T of 17 is moving according to the pseudo-solid strategy. However, the boundary of 
uj is not exactly tracked in either of these meshes. An unstructured mesh consisting of triangles is 
used to approximate 17. The mesh of is obtained by displacing the nodes of T using the discrete 
displacement field Vh which approximates v. 

4.1 Discretization of the state problem 

Let e > 0, -0 e U"'^ be given. To simplify notation, the penalty parameter e at the discrete solution 



will be omitted. The finite element discretization of (Ve{ip) I reads as follows: 
Find {uh,Ph,Vh) G Vhi^^J x Qh 'x Wh such that 

/ Vuh-Viphdx-j (phds+ / G£{'ilj){uh - l)(phdx = Vv?^ G V/i(0„^) 



dn^j 



h 



la{vh):e{wh)dx- lphn-Whds+ lGeiip)vh-Whdx = \/whGWh, 
^ Jn Jdn Jn 

where V/j(ri„^), Wh = Wh x Wh, Qh are finite element approximations of H^{Q.v), [H^{Q)Y and 
L^(90) respectively. Here we shall use linear triangular elements for constructing Vh{^Vh) and Wh, 
while Qh = Whl^fj. 
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4.2 Discrete optimization problem 



In the traditional level set method the function ip is taken to be a function of pseudo-time t, ip := 
ij){x,t), and the optimization process is realized by solving the Hamilton- Jacobi equation 

^+t;„||VV'|| =0. (3) 

Here Vn is the velocity derived by the means of sensitivity analysis, often done on the continuous 
level. The function ijj is then advanced towards the steady state in pseudo-time, see e.g. [101 . 

Despite of the conceptual simplicity it is not so straightforward to implement the conventional 
level set method due to the need of appropriate upwind schemes, an extension of the velocities and 
re-initialization algorithms. Indeed, since the Hamilton- Jacobi equation does not in general admit a 
smooth solution, an appropriate upwind scheme must be used for the time integration. The velocity 
Vn is often meaningful only on the boundary du, and must be extended to the whole domain, or at 
least into a neighborhood of doj. Finally, the function ip should be an approximation of the signed 
distance function, i.e. il^{x) sign(V'(iE))| |a; — a;o||, where xq is the closest point to x for which 
^^{xq) = 0. To force this property a re-initialization procedure is often used. 

Several approaches to overcome these difficulties have been proposed. In lITll radial basis func- 
tions (RBF) are used to define the function ip, and the Hamilton-Jacobi equation is transformed into 
a system of ordinary differential equations. In lfT2l the function i]j is constructed combining param- 
eterized primitives with a radial basis function representation of the so called R-functions. In |(T3]| . 
11.1 4J the function ip is approximated by the same shape functions on the same mesh used to solve the 
state problem. In this paper we follow |[T5l and utilize the compactly supported -continuous radial 
basis functions [ 16| to parameterize explicitly the level set function. Then the geometry will inherit 
smoothness properties of the underlying parameterization. Moreover, we have a fixed set of design 
variables, and one can use sophisticated optimization methods instead of performing integration in 
the pseudo-time. 

We introduce a set of x A'^ basis functions, whose knots are placed in the interior of the domain 
D as follows. The coordinates of the knot are given by 

Xij = Xmin + (i - 1) ^""^ ^^'"^^ ^ i,j = 1, (4) 

, / • -,\ Umax ymin . . ^ ,t 
Vij = Vmin + [l - ^) ]S[ _l ' ^:' = 1'---'^' (5) 

where Xmm and Xj^ax are the minimal and maximal x-coordinates of the rectangle D, respectively 
and similarly ymin and Umax- The RBF associated with this knot is then 

i^ijii^ij) = max{0, (1 - rij)}"^ {Arij + 1), (6) 

where 

I s 

Here > is a given parameter, the radius of the support (see Figure |2]l 
ping of the supports of ipij we define this parameter as 

. / Xyjiax Xmin Umax Umin \ 

's = ^-^^-[ , )■ (8) 



(7) 

To guarantee the overlap- 
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Figure 2: Surface plot of a radial basis function 

Using the assumption that $7 is a star-like domain, the objective functional ([T]l will be given in the 
discrete setting by 

J{cx) := J{vh{cx)) = I \gc.{e) - gt{e)f d9, (10) 
Jo 

where a is the vector of the discrete design variables, V}i{(y.) is apart of the solution to {V^{iI)n{ol))), 
ga{0) is the radius of the free boundary corresponding to a, and gt{0) is the radius of the target 
boundary at the angle 6. The set U"''^ of all admissible level set functions is represented by functions 
of the form ^ with cxeU := [omin, amax]^""^ C M^''^, with Omin, "max given. 
Thus, the finite-dimensional optimization problem to be realized reads as follows: 

Find ct* such that 

(11) 

J{a*)<J{a) yaGU 
subject to {'P^{'tpN{a))). 

4.3 Construction of the initial guess 

Obviously, the location of the free boundary corresponding to a given inner boundary duj is not known 
a priori. Therefore we start the optimization from a simple configuration where the location of the 
free boundary corresponding to the initial uj is known. 
Our initial guess of a is constructed in such a way that 

i^N{xij,yij) ='ipi{xij,yij) = 1,...,N (12) 
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where {xij,yij) are the knots of the radial basis functions and ipi is given by ( |43| ) in Appendix 2 with 
R = 1. We obtain a system of linear equations, which is known to be invertible. The initial values 



for the design vector a are solutions to ( [T2| ). The zero level set of such function ^p then approximates 
dB{0, 1). 

For given 7, the initial reference domain 17 is taken to be the circle B{0, C(l, 7)) (for the defini- 
tion of C(l, 7) see Appendix 2). 



4.4 Algebraic form and Sensitivity analysis 

Let tp E Uad be given. In the pseudo-solid approach we simultaneously seek the scalar function 
u, the press ure p, and the deformation field v which deforms the reference domain O into the one 



that solves (V(uj 



The elasticity system ( P^(^)|) 3 is th us solve d in the undeformed configuration 



Q of the pseudo-solid, whereas equations {V^{ip) ii and (V^{ip) I2 are solved in the deformed one. 



Therefore, they have to be discretized by different meshes, too. Let us denote the nodal co-ordinates 
of the triangulation T of 17 by X. We simply transform this triangulation into the one of 17^^ which 
is characterized by the nodes X = Fh{X), where is defined by the discrete displacement field 
being the approximation of v, i.e. 

X = X + Vh{X). (13) 



(V^{tp))2 resulting from an appropriate discretization can be 
0, respectively and the linear elasticity system (V^{ip) Is as 



The algebraic form of {V^{ip) 

written as ri(q„, g^) = 0, r2(q„, q„) 

rs{q^, Qp) = 0, where and Qp are the nodal values of Uh, and ph, respectively. Here the 

dependence of ri and r2 on is through the nodal co-ordinates, as specified by ( [T3] ). Dimensions 
of the vectors q„, q^ and q^ are n, 2n and rig respectively, where n is the number of the nodes in T 

and rie is the number of the nodes on 917. 
Let us introduce the following notation: 



Qv 

L<lp. 



and 



ri 



(14) 



Then the algebraic form of the discretized coupled system (V^{ip) I can be written in short as r{q) 
0. This system will be solved using Newton's method: 



/ gr(gW) Y 



(15) 



with the Jacobian matrix 



dr 
dq 



dri 

din 

dr2 





dr2 






dqp. 



(16) 



Remark 1 The discrete cost function J' defined by ([TO]) is not differentiable due to the discontin- 
uous Heaviside function H. If one wishes to use descent type optimization methods (and this will 
be our case), then smoothing of the Heaviside function is necessary (see Appendix 1) to make the 
components of the residual vector r continuously differentiable functions of the design variables. 
Thus, while assembling the discrete system arising from (V^{Tp) ), H will be replaced by a smoothed 
function Hp, i.e. Ge = -Hp in (V^{ip) I in what follows. 
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Remark 2 Notice that the equations are coupled in a quite comphcated way. For example, the resid- 
ual ri depends naturally on q„, but also on q„ through the shapes of the elements as specified by ( [T3| ). 
Moreover, since the mesh is moving and ^p is a function of location, ri depends on also through 
Hji^tp). This dependence must be taken into account especially during the sensitivity analysis phase 
in order to obtain perfectly consistent derivatives. 

Despite this nonstandard coupling between the equations, the Jacobian matrix ([16]) is easy to 
compute using the sparse forward mode automatic differentiation ifTTl . Our implementation of the 
automatic differentiation technique is described in ifTSl . For a general introduction to the principles 
of automatic differentiation see |fT9l . 

If J is smooth then using the well-known adjoint approach, the gradient VaJ can be computed 
from 



V.J =13^1 (17) 



where the adjoint vector u solves the adjoint equation 



dr \ 

\ V=-VqJ. (18) 



dq) 



The required Jacobian matrices and gradients in ( 17 1 and ([18]) can be again easily computed using the 



tools of automatic differentiation. Notice that the Jacobian dr/da has a sparse structure, since the 
radial basis functions are compactly supported. This sparsity is automatically exploited, since we use 
the sparse forward mode automatic differentiation technique. To avoid going through all radial basis 
functions while evaluating ^, a quadtree data structure is exploited to exclude RBFs that can not have 
a non-zero value at the point of evaluation. 

4.5 Optimization strategy using remeshing 

As explained in the previous paper |[T], remeshing is sometimes needed since the mesh deformation 
approach can not handle too large displacements. Indeed, the Newton method used to solve the 
coupled system may not converge in case of excessive deformations of Q. Moreover, if the mesh gets 
too distorted, significant errors in the solution and numerical instabilities may appear. 

In this paper we adopt the following adaptive optimization strategy: to avoid too large deformation 



fields Vh in (T'^(V') I the reference domain is re-initialized after /cmax optimization steps. Moreover, 
the half width 5 of the gray region related to the smoothing of the Heaviside function (see Appendix 1) 
is determined adaptively. This 5 is then used to determine the smoothing parameter /3, i.e. /3 := I3{5). 
Note that due to this adaptive strategy /3 varies spatially, too. Therefore there are two smoothing 
parameters e, 5 related to the discrete pseudo-solid problem which we thus denote by {Vj^^^t/j)). 

The strategy of choosing 6 is also altered after A;max optimization steps. In early steps a larger 
value of 6 is used, which makes the objective function smoother and enables fast progress. The value 
of 6 is then gradually decreased, improving the approximation of the exact Heaviside step function. 
We choose the half width 6 of the desired gray region after i:th re-initialization step to be 

6 = liim.. - i + 1)K (19) 

where imax is the number of re-initializations to be done and h is the characteristic mesh size. 

We choose the number of re-initializations to be modest, e.g. imax = 8. In the early stages of 
optimization the progress is rapid, and changes in the design domain u are large. Therefore we let the 
number of steps between re-initializations to increase in the course of optimization, e.g. k^ax = 5-2*. 
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The re-initialized optimization process can be described as the following algorithm: 



Algorithm 1 {imi,^,h,toler given) 

Compute initial guess a using eq. (12) ; (20) 

SetJ7 := 5(0,C7(1,7)) ; (21) 

fori = l,...,w (22) 

Generate triangulation T of 17 ; (23) 

Set 6 := ^(Wax -i + l)h and /c^ax = 5 • 2* ; (24) 

do (25) 

Set k:=l; (26) 

Solve state problem {V^^^^ipNia))) ; (27) 

Evaluate J{a) and VcJ{a) using (fTOjl, ([17) and ([T8]l; (28) 

Find descent direction a of ^ at a ; (29) 

Update design a := ot + a ; (30) 

Set A: := A; + 1 ; (31) 

while k < /cmax and ||q!|| > toler ; (32) 

Reinitialize J7 := using the latest Vh ; (33) 

end (34) 



The re-initialization of U in step p3] ) is done as follows. One fits by least squares a cubic B- 
spline curve to the outer boundary of the deformed mesh T{vh) obtained as the solution to the free 
boundary problem corresponding to the current design a. Nodes are then distributed on this curve, 
and this geometry is given to the mesh generator Gmsh ||20| . 

Thus, in algorithm (pT])-((34l) we basically solve i^ax different optimization problems, performing 
a mesh regeneration for each problem. The initial guess for the next problem is always the best design 
found so far, i.e. the end result of the previous problem. 



5 Numerical examples 

In this section we illustrate the performance of the proposed method. The target domains are the same 
as in im. The value of the penalty parameter e is 10~^. We used the gradient based optimizer Donlp2 
ETll to realize the step (|29]l in Algorithm 1. The parameters amin, ctmax were chosen as ±10^'', i.e. 
we have practically unconstrained minimization problem. The performed numerical computations 
indicated no need to pose more strict constraints. 

Example 1 The target Ft is the "rounded square". The length of sides of the square is 4. Each corner 
is rounded using a quarter of a circle of radius 1 . For the magnitude of the normal derivative the value 
7 = — 1 was used. It was observed in HI that if the family O of admissible inner inclusions contains 
only star-like domains uj, then the target Tt can never be matched. Moreover, if the number of the 
design variables increased, the boundary dio became more and more oscillating. 

The final reference domain Q was discretized using 18842 linear triangular elements. The number 
of RBFs used was 30 x 30. The initial value of the cost functional was 1.09. After 630 optimization 
steps (and 2376 function evaluations) the cost function reduced to the value 3.78 • 10~^. 
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The final geometry and the contour plot of His{ip'^) are shown in Figure [3] 




Figure 3: Left: the final geometry and Hp{il)*j^). Right: zoom of the framed area. 

In addition to "black" (i.e. w) and "white" (i.e. regions, smoothing of the Heaviside function 
also produces by its construction "grey" regions (where < Hp{ip*j^) < 1) near the zero level set of 
representing duj. In the final design however grey regions which can not be interpreted as being 
close to du) may appear. In Figure |3] we see one of those problematic regions. By examining the 
potential (see Figure [4]) in that area, we find that the maximum value of Uh in this region is only 
about 0.86, whereas in lo we should approximately meet the condition = I- 




Figure 4: Contour plot of the potential Uh related to Example [T] 

Example 2 To prevent such problematic grey regions we propose to add the following penalty term 
to the cost functional J': 

Jr, = V [ H2^W{u - if dx, 7? > 0. (35) 

Jo. 
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The idea behind is the following. The penalty Jr^ will be zero where H213 is zero or u = 1 is met. 
Near the boundary of a; the solution u is close to 1 , and J^^ is small. However, in regions like the one 
shown in Figure [3] the value of is larger since Uh is much less than 1 . Notice that instead of the 
original value used in the state problem, a two times larger value 2/3 was used in the smoothing of the 



Heaviside function in (35 1 



This approach turns out to be effective. In Figures [5] and |6] we show the results for the same 
problem as in Example [T] but with the penalty term included. The value of the penalty parameter 
was r/ = 0.1. Now the problematic grey regions near the corners have disappeared, and have been 
replaced by regions where the Heaviside function has the value Hp = 1. Also the potential is 
close to 1 near this region (see Figure |6]). This time the optimizer needed 1009 iterations and 4098 
function evaluations. The final values are = 7.67 - 10^^ and J„ = 8.42 • 10^^. 




Figure 5: The final geometry and Hi3{ip^). 




Figure 6: The potential Uh related to Example [2j 

Example 3 Let 7 = — 2 and let the radius function defining the target boundary Tt be chosen as 
follows: 

g^(e) = 0.5 cos(6') + 0.8 003(26*) + 2, [0, 27r[. 
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The final finite element mesh consists of 15820 linear triangular elements, and the number of RBFs 
is 15 X 15. Value of the penalty parameter was r] = 0.01. The initial value of the objective function 
(i.e. the sum of ^ and ^ri) was 4.77. After the total of 356 optimization steps and 1559 function 
evaluations this value was reduced to 1.90-10"^. The results of computations are depicted in Figures|7] 
and [8] 

This problem was solved in |^ using a fixed topology approach. The results of computations 
led to a conclusion that the inner boundary consists of more than one component. From this rea- 
son two holes as an initial approximation of the inner boundary were introduced "by hand", each 
parameterized by the radial co-ordinates. The obtained result was similar to that in Figure [7] 




Figure 7: The final geometry and Hp{^l)*j^) 




Figure 8: The potential Uh related to Example [3j 
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6 Conclusions 



In this paper we have considered topology optimization with the state constraint given by a free 
boundary problem of Bernoulli type. To solve efficiently the free boundary problems during the 
optimization, the pseudo-solid approach is applied. Its main advantage is that there is no explicit 
parametrization of the shape of the free boundary using e.g. splines. The novelty of the numerical 
method proposed in this paper is the combination of the pseudo-solid approach to tackle the free 
boundary problem with a parameterized level set method for shape optimization. It has been found 
already in [1] that the problem is very badly conditioned as many different choices of w may lead 
to nearly identical free boundaries. Therefore the progress of the optimization is often slow. The 
proposed method can be applied in an analogous way to topology optimization problems governed 
by other free boundary problems. 

References 

[1] J. I. Toivanen, J. Haslinger, and R. A. E. Makinen, "Shape optimization of systems governed by 
Bernoulli free boundary problems," Computer Methods in Applied Mechanics and Engineering, 
vol. 197, pp. 3803-3815, 2008. 

[2] M. Flucher and M. Rumpf, "Bernoulli's free boundary problem, qualitative theory and numerical 
approximation," J. Reine Angew. Math., vol. 486, pp. 165-204, 1997. 

[3] K. Karkkainen and T. Tiihonen, "Free surfaces: shape sensitivity analysis and numerical meth- 
ods," International Journal for Numerical Methods in Engineering, vol. 44, no. 8, pp. 1079- 
1098, 1999. 

[4] C. Cuvelier and R. M. S. M. Schulkes, "Some numerical methods for the computation of cap- 
illary free boundaries governed by the Navier-Stokes equations," SIAM Review, vol. 32, no. 3, 
pp. 355-423, 1990. 

[5] J. Haslinger and R. A. E. Makinen, Introduction to Shape Optimization: Theory, Approximation, 
and Computation. SIAM, Philadelphia, 2003. 

[6] R. A. Caimcross, P. R. Schunk, T. A. Baer, R. R. Rao, and P. A. Sackinger, "A finite element 
method for free surface flows of incompressible fluids in three dimensions. Part I. Boundary 
fitted mesh motion," International Journal for Numerical Methods in Fluids, vol. 33, pp. 375- 
403, 2000. 

[7] M. Souli and J. P. Zolesio, "Arbitrary Lagrangian-Eulerian and free surface methods in fluid 
mechanics," Computer Methods in Applied Mechanics and Engineering, vol. 191, pp. 451^66, 
2001. 

[8] S. Osher and J. A. Sethian, "Fronts propagating with curvature dependent speed: Algorithms 
based on Hamilton- Jacobi formulations," Journal of Computational Physics, vol. 79, pp. 12-49, 
1988. 

[9] S. Osher and R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces. Springer Verlag, 
2003. 

[10] G. Allaire, F. Jouve, and A.-M. Toader, "Structural optimization using sensitivity analysis and 
level-set methods," Journal of Computational Physics, vol. 194, pp. 363-393, 2004. 



14 



[11] S. Wang and M. Y. Wang, "Radial basis functions and level set method for structural topology 
optimization," International Journal for Numerical Methods in Engineering, vol. 65, pp. 2060- 
2090, 2006. 

[12] J. Chen, V. Shapiro, K. Suresh, and I. Tsukanov, "Shape optimization with topological changes 
and parametric control," International Journal for Numerical Methods in Engineering, vol. 71, 
no. 3, pp. 313-346, 2007. 

[13] T. Belytschko, S. P. Xiao, and C. Parimi, "Topology optimization with implicit functions and 
regularization," International Journal for Numerical Methods in Engineering, vol. 57, pp. 1 177- 
1196, 2003. 

[14] P. Neittaanmald, A. Pennanen, and D. Tiba, "Fixed domain approaches in shape optimization 
problems with Dirichlet boundary conditions," Inverse Problems, vol. 25, no. 5, 2009. Article 
number 055003 (18pp). 

[15] Z. Luo, M. Y. Wang, S. Wang, and P. Wei, "A level set-based parameterization method for 
structural shape and topology optimization," International Journal for Numerical Methods in 
Engineering, vol. 76, pp. 1-26, 2008. DOT: 10.1002/nme.2092. 

[16] H. Wendland, "Piecewise polynomial, positive definite and compactly supported radial functions 
of minimal AegvQQ" Advances in Computational Mathematics, vol. 4, no. 1, pp. 389-396, 2005. 

[17] C. H. Bischof, P. M. Khademi, A. Bouaricha, and A. Carle, "Efficient computation of gradients 
and Jacobians by dynamic exploitation of sparsity in automatic differentiation," Optimization 
Methods and Software, vol. 7, pp. 1-39, 1996. 

[18] J. 1. Toivanen and R. A. E. Makinen, "Implementation of sparse forward mode automatic dif- 
ferentiation with application to electromagnetic shape optimization," Optimization Methods and 
Software, 2010. To appear, published online. 

[19] A. Griewank and A. Walther, Evaluating Derivatives: Principles and Techniques of Algorithmic 
Differentiation. SIAM, Philadelphia, 2nd ed., 2008. 

[20] C. Geuzaine and J.-F. Remade, "Gmsh: A 3-D finite element mesh generator with built-in pre- 
and post-processing facilities," International Journal for Numerical Methods in Engineering, 
vol. 79, pp. 1309-1331,2009. 

[21] P. Spellucci, "An SQP method for general nonlinear programs using only equality constrained 
subproblems," Mathematical Programming, vol. 82, pp. 413^148, 1998. Software available at 
http://plato.la.asu.edu/donlp2.html. 

Appendix 1 

Smoothing the Heaviside function 

The expression for the C^-smoothed Heaviside function is done by 



0, 
1, 



ify < -/3 
ify >/3 



(36) 




otherwise. 
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where /? > is a given constant. 
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Figure 9: The graph of H/^, (3 



lfip{x) /3[ then a; belongs to the so-called grey region, i.e. Hp{4'{x)) g]0, 1[. Such regions 

appear near the boundaries of the design domain uj. To maintain a good quality of the gradients 
that we obtain, we would like to have the geometrical width of the grey region to be approximately 
constant everywhere. If the scalar function ip approximated the signed distance function, this would 
automatically be the case. However, we do not make this requirement, but instead we propose the 
following procedure. 

Let 25 be the desired width of the grey region. Consider a point on the zero level set, and a 
point X located ntx = xq-\- 6n, where n is the unit normal vector to the zero level set of tp. The unit 
normal vector is n = VV'/||VV'||, thus we can write 



X = xq + 5n = xq + 6 



Vijixo) 



\\Vij{xo)W 

Since on the other hand we have the first order approximation 

V^(a;) ^ tp{xo) + {x- xo) ■ Vij{xo), 



we get from ( 37 1 and ( 38 1 that 



Thus if we use 



\\V^{xo)\\ 



f3{x) = 6\\V^{x)\\ + 10- 



(37) 



(38) 



(39) 



(40) 



as the Heaviside parameter, the width of the grey region will be approximately 25 everywhere. The 
constant 10^® is added to prevent division by zero, since /3 appears in the denominator in ([36]l. 
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Appendix 2 

Analytical solution in the circular domain 

Let uj = B{0, R), where B{0, R) is the circle of radius R centered at the origin, and C := C{R, 7) 
be a constant such that 

C ln(C) - C ln(R) = - i , 7 < 0. (41) 

7 

The function 

f{R, 7) := C7 ln( + y2) _ c-f ln{R) + 1 (42) 

satisfies A/ = in J] \ cJ, / = 1 on Scj and V/ • n = 7 on dB{0, C{R, 7)). Thus, giv en 7 < 
and CO = B{0,R), we know that the analytical solution of the free boundary problem {V{lo)) is 



u = f{R, 7) and $7 = -6(0, C{R, 7)). Finally, a level set function representing B{0, R) is given by 

4^r{x, y) = R- ^x'^ + y2. (43) 
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